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BAYESIAN PREDICTION OF MEAN SQUARE ERRORS WITH 

COVARIATES 



by P. A. Jacobs and D. P. Gaver 

1. INTRODUCTION AND SUMMARY 

Numerical meteorological models are used to assist in the prediction of 
weather. Each run of a numerical model produces forecasts of meteorological 
variables which are used as preliminary predictions of future values of these 
variables. These initial predictions are referred to as first-guess values. In 
this paper first-guess values will refer to the most recent 12-hour forecasts. 

In certain areas of the world, observations of forecasted variables become 
available. Prior to the next run of the numerical model a multivariate 
optimal interpolation analysis updates a first-guess value of a variable by 
adding to it a weighted observed value of the variable if it is available. The 
weight multiplying the observed value depends on estimates of the mean 
squared error of the first-guess value and the mean squared error of the 
observation; cf. Goerss et al., [1991, a, b]. Thus it is of importance to predict 
such first-guess squared errors. 

The general problem of modeling and predicting mean square errors is 
important but not widely studied; see Davidian and Carroll (1987), Nelder and 
Lee (1992), Aitken (1987), McCullagh and Nelder (1983). 

In Jacobs and Gaver (1991, 1992) statistical models for the error of the first- 
guess are used to predict mean square error for first-guess wind components. 
The models assume that the error of the first-guess is normal with mean 0 



1 



and variance which is a function that is log-linear with suitable covariates. 
The cross-validation results of those papers suggest that covariates do have 
some predictive ability for the mean square errors. However, the relations 
change over time. 

In this paper we introduce a procedure for recursively updating the 
estimated parameters of the variance function. The approach is Bayesian 
with recursive updating using an approximation based on the Laplace 
method; cf. deBruijn (1958). 

In the next section the model is introduced. Details of the updating 
procedure are also given. 

The third section presents results of using the procedure to predict mean 
square first-guess wind component errors. The data consist of measurement 
and 12 hour forecasts (first-guess values) of u and v wind components at 850 
mb, 500 mb, and 250 mb pressure levels from 93 stations in North America, 
25N-75N. The measurement values (if available) are subtracted from the 
first-guess values to obtain observations of the first guess error. The 
covariates considered are wind speed and resultant wind, (the sum of the 
squared difference of the w— wind component at two consecutive 12 hour 
periods and the squared difference of the t;-wind component at the same 
times). The resultant wind is a measure of the change in the atmosphere. 
Higher wind speeds suggest more activity in the atmosphere. 

The results of the data analysis suggest that the covariates do have 
predictive ability. The models using observed wind speed and resultant wind 
have more predictive ability than those using the first-guess values of wind 
speed and resultant wind. Further, models that use both wind speed and 
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resultant wind have more predictive ability than those using either one by 
itself. The change of the model parameters with time appears to be slow. 
This suggests that while the relationship of the mean square error and the 
resultant wind and wind speed is changing, it may not be necessary to update 
the model parameters in every period. 

2. THE MODEL AND STATISTICAL PROCEDURE 

Let Yiit) denote the first-guess error at location i at time t; i = 1, ..., L. Let 
Xi(f) = = 1, ..., p) denote the covariates at location i at time f. 

Consider the following model for the first-guess errors. 

p{}^(r)ecfy,.(r);i = l,...,L|x,(r),[3(r),y/4x^(45<t,7 = l,...,L} ( 2 . 1 ) 

where x,(f)P(f) = ^^■^^Xij{t)Pj{t);ihai is, given P(f), the first-guess errors are 
conditionally independent normally distributed with mean 0 and log linear 
variances 

( 7 ?(f) = exp{x,-(t)p(/)} (2.2) 

independent of everything else. 

The coefficients P(f) are modeled as changing according to a random walk 

P(f-H) = p(f)+a>(f + l) (2.3) 

where {co(f)} are independent multivariate normal random variables with 
variance-covariance matrices {W(f)). The matrix W(t+1) is independent of 
{y,(s), x,(s), p(s), s<t I = 1, ..., L} 
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In the next subsection we suggest a Kalman filter-like procedure to 
produce successive estimates of P(f) as new data become available. The 
procedure is based on a Laplace approximation to an integral. 

2,1 An approximate Updating Procedure 

Assume the posterior distribution of P(f) given (y/fs), / = 1, ..., L, s < f} is 
multivariate normal with mean m(0 and variance-covariance matrix E(f). 
Since it is known that 



the prior distribution of P(f+1) is multivariate normal with mean m(t) and 
variance-covariance matrix 



A description of the procedure used to determine W(f+1) appears in Section 3. 

The forecast/ prediction distribution of {Y,(f+1); i = 1, L) in terms of data 
up to time t and covariate values at time t +1 is 



P(f+1) — P(0 + co(f+l). 



R(f) = Z(0 + W(f+1) . 



(2.4) 




(2.5) 



where 



/p(f+l)(b) = ((2^f detR(f)) 2exp|-|(b-m(f))R \f)(b-m(f)) j. (2.6) 
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We now approximate the integral by the Laplace method; cf. Easton 
(1991), Cox and Hinkley (1974), de Bruijn (1958). Let the exponent of the 
integrand be 




L 



+ yf exp{-x,b} + (b - m(f))R ^ {t)(b - m(0)' 

. 1=1 



+ K 



(2.7) 



where X is a constant and we let Xf = x,(t+l). Differentiating, we obtain 




1 - y} exp{-x,b}] + 2R^-.^ (t)(b - m{t)) 



(2.8) 



where Rj^it) denotes the row of R \f) and 



a" ..__i 

dbjdbk ^ 2 



L 

Xy? exp{-x,b}xyx,-;t + 2RJ^^ 

Li=l 



(2.9) 



Use a Newton procedure to solve the system of equations 






/ = l,...,p- 



(2.10) 



for m(f+l). Solve for X(f+1) using the second derivatives of g evaluated at 
m(f+l); that is X(f+1) is minus the inverse of the matrix whose (/, k) entry is 

The posterior distribution of p(f+l) given {YKs), x,(s), i = 1, ..., L, s = 1, ..., f+1} is 
approximated by a multivariate normal distribution with mean m(f+l) and 
variance X(f+1). The estimate of P(f+1) is m(t+l). 

The predicted mean square error for the ith location at time f+2 is 
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exp^ 






l/=l 



3. DATA ANALYSIS 

In this subsection we report results concerning using regression-like 
models for the mean square error of the first-guess with recursively updated 
parameter estimates to predict future mean-square errors of the first-guess. 

The data consist of measurement and 12 hour forecasts (first-guess 
values) of u and v wind components at the 850 mb, 500 mb and 250 mb 
pressure levels from 93 stations in North America 25N-75N for the month of 
July 1991. The forecasts are produced using the NOGAPS Spectral Forecast 
Model; cf. Hogan et al., (1991). Each station has measurement and first-guess 
values for every 12 hours; there are some missing observations and 
suspicious values of wind components equal to 0. These missing and 
questionable values are deleted from the data set. The measurement values 
(if available) are subtracted from the first-guess values to obtain observations 
of the error of the first-guess value. 

Let U(0;t), (respectively P(0;f)), be the observed w-wind, (respectively 
y-wind) component at time f. Let (respectively V(f;t)), be the first-guess 

u-wind (respectively i^-wind) component at time t; Uif;t) is the forecasted 
value of U(t) made 12 hours previously. The first-guess error for the w-wind 
component is 

Yit) = Uif;t) - LK0;f). (3.1) 

The following covariates are considered in the log-linear model for the mean 
square error of the first-guess. 
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r(0; t) = [(a(0; t)-U{0;t- if + [V{0;t) - V{0; t - if 



(3.2) 



s{0;t) = \u{0;f+V{0;ty 



(3.3) 



K/; 0 = (^(/; 0 - ^(/; t - i)f + - v{f; t - i)f 



(3.4) 



s{f;t) = [u{f;f^V{f;tf 



1 _ 

2 



(3.5) 



The resultant wind r(0;t), (respectively r(f;t)), is a measure of the observed 
(respectively forecasted), change in the wind. The variable s(0;0, (respectively 
s(/;/)), is the observed, (respectively forecasted), wind speed. Higher wind 
speeds suggest more activity in the atmosphere. 

3.1 The Models 

The following models for the mean square error are considered 
One Variable Models: Observed Covariates 

1. Given Poit), P\(t), the first-guess errors at each location {y,(0; i = 1,. ., L) 
are independent normally distributed with mean 0. The variance of 
Yi(t) is the following function of the observed resultant wind at 
location i at time t, r,(0;/) 

cr|(l;f;r,•fO;^;) = exp{/?o(0 + /?i(Or,(0;O} (3.6) 

where rj(0;f) is the observed resultant wind at location i at time t and 
Yj(t) is the first-guess error at location i at time t. 

2. Given Pq{ 1), P\{t), [Y jit); i = 1, ..., L} are independent normally 
distributed with mean 0. The variance of Yi(t) is the following 
function of the observed wind speed at location i at time t, s,(0;f) 
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(3.7) 



a 2 (l;t;si( 0 ;t)) = exp{)3o(0 + )Si(Os,(0;f)} 
where s,(0;f) is the observed wind speed at location i at time t. 
Two-variable Model, Observed Covariates 

3. Given j3o(0, PiU), (TjCf); i = 1, L} are independent normally 
distributed with mean 0. The variance of Y,(f) is the following function 
of both the resultant wind and wind speed at location i at time t 

C^it; n(0;t); si(0;t)) = exp{/3o(f) + j3i(f)rf(0;f) + p2(Os,(0;f)}. (3.8) 

Similar one-variable and two-variable models but using first-guess values of 
the covariates are also considered. In all cases the first-guess error and the 
covariates are all evaluated at the same pressure level. 

The regression parameters P(f) are assumed to evolve according to the 
random walk given by (2.3) 

3.2 The Data Analysis 

The results of Jacobs and Gaver (1992) suggest that of the models using 
observed values for covariates, the models for the 850 mb pressure level have 
the most predictive value. It is also suggested that of the models using first- 
guess values for covariates, the models for the 250 mb pressure level have the 
most predictive ability. As a result in what follows we will restrict our 
attention to these two cases. 

a. Estimation and Prediction of 850 mb First-guess Mean Square Errors using 
Observed Wind Covariates. 

The estimation procedure described in Section 2 was used to recursively 
estimate the regression parameters P(f) for each of models (3.6) - (3.8) for 850 
mb first-guess errors using observed 850 mb wind covariates. The initializing 
estimates of P(f) are the estimates obtained using all April data recorded in 
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Jacobs and Gaver (1992); the initial variance-covariance matrix is the identity 
matrix. The estimates from April are used since April appeared to have more 
predictive ability for July than February, cf. Jacobs and Gaver (1992). 

The variance-covariance matrix of the innovation Wt of the random 
walk (2.3) was taken to be a constant times the identity matrix. Preliminary 
explorations based on values of the predictive log-likelihood using different 
values of the constant suggest that for purposes of prediction, the constant 
should be very small. In what follows we set the constant equal to 0. 

Figure 1, (respectively Figure 2) presents plots of the estimates of the 
slopes, e.g., p-[(t) and ^ 2 ( 0 , as a function of time which is labeled 1, 2, .... for u- 
wind (respectively n-wind) component error. 

Figure 1 presents the values of estimates of the parameter multiplying 
r,(0;f) (respectively s,(0;t)) for the one parameter models (3.6) and (3.7) in the 
upper graph. The lower graph presents the values of the estimates of the 
parameter multiplying r,(0;t), (respectively s,(0;t)) for the two-variate model. 
There appears to be a slight trend in the estimates. 

Figure 2 presents the values of the estimate of the parameters multiplying 
ri(0;t), (respectively Si(0;t)) as o, (respectively +), for the respective one-variate 
models. The values of the estimates of the parameters multiplying r,(0;f), 
(respectively Si(0;t)) for the two- variate model are presented as x (respectively 
V) for each time t. These graphs suggest more evidence of a trend in the 
estimates. Note that the estimates of the slopes are positive. Hence increased 
values of the resultant wind r,(0;t) and/or wind speed s,(0;f) are associated 
with increased variance of the first-guess value. This is plausible physically. 
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since a large value of ri(0;t) is indicative of a change in the atmosphere and a 
large value of Si(0;t) is indicative of greater activity in the atmosphere. 

To assess the predictive ability of the models, the models with parameters 
estimated at time t are used to forecast the variances of the first-guess errors at 
time f+1. 

One procedure to informally assess the predictive ability of the models is 
by binning the data. To assess models (3.6) and (3.8) the data (yj(0, r,(0;f), 
Si(0;O) are binned into 10 bins based on ordering the values of r,(0;f) for all 
time t from smallest to largest. The data in the first bin correspond to the 
smallest values of r,(0;0; the data in the 10th bin correspond to the largest 
values of r,(0;0- Each bin contains about 1/lOth of the data with the 10th bin 
containing a few more data. The averages of the predictive variances for 
models (3.6) and (3.8) are computed for each bin. The average is also 
computed for each bin. 

To assess models (3.7) and (3.8) the same procedure is used but the 
binning is based on the values of s,(0;f). 

Figures 3 and 4 present graphs of the log[average y,(f)^] in each bin versus 
log [average predictive variance] in each bin for models (3.6) and (3.8) and 
models (3.7) and (3.8). If a model were perfect, the points should be close to 
the 45° line shown. 

The figures suggest that of the two one-variate models, the one using the 
resultant wind r,(0;f) has the better predictive ability. The two-variate model 
appears to have similar predictive ability to the model using only r,(0;0. 

Another procedure to assess predictive ability is to compute the log- 
likelihood using estimated values of P(f) (eg. m(0) in the term for the first 
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guess errors at time f+1. Larger values of the (predictive) log-likelihood 
indicate better predictive ability. The (predictive) log-likelihood up to 
addition and multiplication by constants is 

T-\ 

'=-1 
( = 1 

Table 1 presents values of / for the one-time step ahead predictions of the 
variance; these values appear in the column Iterative. Also displayed are the 
values of / obtained by estimating the parameters once using all the data; this 
value of / is a goodness-of-fit value and appears in the column labeled All; 
the estimates used to obtain the goodness-of-fit value of / are those appearing 
for July in Table 7 of Jacobs and Gaver (1992). 

Four models are considered: constant variance (no dependence on 

variables), two one-variate models (3.6) and (3.7) and the two-variate model 
(3.8). The parameter of the constant variance model using all the data is 
estimated using maximum likelihood; this estimate is used to calculate the 
goodness-of-fit value of / for the constant variance model. 

The goodness-of-fit constant variance value of I is smaller than the 
prediction values of / using models with covariates. This behavior indicates 
that the covariates do have predictive ability. The prediction value of / for 
the two-variate prediction model is larger than that for either one-variate 
model indicating that both covariates have some predictive ability. Of the 
two one-variate models, the one using ri(0;t) has the larger prediction value 
of /. The closeness of the prediction values of I obtained by iteratively 
updating the estimates and using them to predict variance of the next time 
period and the goodness-of-fit values of I obtained by using model 



^ X, (r -I- l)m(r) + (r -I- l)exp{-x, (/ + l)m(r)} 






L<=i 



1=1 
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parameters estimated from all the data suggest that the updating procedure is 
doing very well. Note that for the one-variate model using ri(0;t), the 
prediction values of / are larger than the goodness-of-fit values of I ; this 
suggests that there is systematic change in the values of p over time for this 
model. 



TABLE 1 

VALUE OF LOG-LIKELIHOOD OBSERVED COVARIATES, 850 mb 



Data 


Model 


Iterative 

(Prediction) 


All 

(Goodness-of-Fit) 


u-wind 


Constant Variance 
One-variate 


-9610.7 


-9510.0 




r(0;0 


-8544.1 


-9098.3 




s(0;f) 


-8607.9 


-8596.2 




Two-variate 


-8489.2 


-8467.5 


y-wind 


Constant Variance 
One-variate 


-9692.7 


-9620.2 




r(0;f) 


-8730.9 


-9294.1 




si0;t) 


-8769.2 


-8756.5 




Two-variate 


-8663.8 


-8652.4 



b. Estimation and Prediction of First-guess Mean Square Errors using First- 

guess Wind Co variates 

The recursive estimation procedure in Section 2 was used to estimate the 
regression parameters P(f) for each of models (3.6) - (3.8) for 250 mb first-guess 
errors using first-guess 250 mb wind covariates; that is, the first-guess wind 
speed at location i at time t at the 250 mb level, Sj(/;f) replaces Si(0;t), etc. The 
initializing estimates of P(f) are the estimates obtained using all April data 
recorded in Jacobs and Gaver (1992); the initial variance-covariance matrix is 
the identity matrix. The estimates from April are used since April appears to 
have somewhat more predictive ability than February for July; Jacobs and 
Gaver (1992). 
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Once again, preliminary exploratory work using the resulting value of the 
predictive log-likelihood indicates that setting the variance-covariance matrix 
of the innovation of the random walk equal to 0 gives the best predictions. 
This suggests that the change in the relationship of the mean square error and 
the covariates is slow. 

Figure 4, (respectively Figure 5) presents plots of the values of the 
estimates for the 250 mb u-wind component errors and the 250 mb z;-wind 
component errors. The values of the estimates multiplying rif;t), 
(respectively s(f;t)) are represented by o, (respectively +) for the one-variate 
models. For the two-variate model, the estimates multiplying r(/;f), 
(respectively s(f;t)) are presented as x, (respectively V). The figures suggest 
evidence of a trend in the estimates. Note that once again all the estimates 
are positive. Thus, increased first-guess resultant wind and/or wind speed 
tends to increase the mean square error 

To assess the predictive ability of the models, the models with parameters 
estimated at time t are used to forecast the variances of the first-guess errors at 
time f+1. 

One procedure to informally assess the predictive ability of the models is 
by binning the data. The data are binned as in the previous subsection. 
Figures 7 and 8 present graphs of the log [average yi(f)^] in each bin versus log 
[average predicted variance] in each bin for models (3.6) and (3.8) and models 
(3.7) and (3.8). If a model were perfect, the points should be close to the 45° 
line shown. 

Table 2 presents values of /, given by (3.9), for the one-time step ahead 
prediction of variance; these values appear in the column iterative. Also 
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displayed are the values of / obtained by estimating the model parameters 
once using all the data; this value of / is a measure of goodness-of-fit and 
appears in the column labeled All; the estimates used in regressions with 
covariates for the goodness-of-fit evaluation of / are those appearing for July 
in Table 11 of Jacobs and Gaver (1992). The constant variance estimate using 
all the data is the maximum likelihood estimate. The results for four models 
are presented; constant variance models, one-variate models (3.6) and (3.7), 
and the two-variate model (3.8). 

Note that all of the iterative prediction values of / for the models with 
covariates are larger than those for the constant variance goodness-of-fit 
value; this suggests that the covariates have some predictive value. 

The iterative prediction value of / for the one-variate model using r(/;f) is 
greater than the goodness-of-fit value obtained by using parameters estimated 
using all the data; this suggests that there is a systematic change in the model 
parameter values over time. For the other regressions using covariates the 
prediction values of / are smaller than their corresponding goodness-of-fit 
values but not by much. The model that maximizes the prediction values of 
/ is the two-variate model suggesting that both covariates have some 
predictive value.. 

Figures 9 and 10 and Table 3 present results for models of the variance of 
the first-guess error at 500 mb level using 500 mb first-guess covariates. The 
implications of the results are similar to those of the 250 mb results. The 
predictive ability of the recursively estimated two-variate model appears 
somewhat greater at the 250 mb level than the 500 mb level; this conclusion is 
based on the values of (/2 - lc)/^c where 1 2 is the prediction value of I for the 
two-variate model and Ic is the goodness-of-fit value of I for the constant 
variance model; the value of this fraction is larger for 250 mb than for 500 mb. 
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TABLE 2 

VALUE OF LOG-LIKELIHOOD FIRST-GUESS COVARIATES, 250 mb 



Data 


Model 


Iterative 

(Prediction) 


All 

(Goodness-of-Fit) 


u-wind 


Constant Variance 
One-variate 


-14,589.4 


-14,516.5 




rif;t) 


-13,640.6 


-14,474.7 




s(/;0 


-13,598.3 


-13,575.5 




Two-variate 


-13,599.5 


-13,579.4 


y-wind 


Constant Variance 
One-variate 


-14,429.1 


-14,363.1 




r(/;f) 


-13,389.2 


-14,270.1 




s(/;f) 


-13,388.4 


-13,363.8 




Two-variate 


-13,351.1 


-13,329.6 



TABLE 3 

VALUE OF LOG-LIKELIHOOD FIRST-GUESS COVARIATES, 500 mb 



Data 


Model 


Iterative 

(Prediction) 


All 

(Goodness-of-Fit) 


u-wind 


Constant Variance 
One-variate 


-11,237.3 


-11,213.5 




r{f-,t) 


-10,561.2 


-11,176.8 




s(/;f) 


-10,518.7 


-10,499.6 




Two-variate 


-10,520.0 


-10,489.0 


i;-wind 


Constant Variance 
One-variate 


-11,204.0 


-11,160.9 




r(f;t) 


-10,431.0 


-11,118.2 




s(/;f) 


-10,423.0 


-10,411.9 




Two-variate 


-10,411.8 


-10,393.3 
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Figure 1 



850 MB V W!ND;RECURSIVE ESTIMATES;OBS WINDS 
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Figure 2 



850 MB U WIND;RECURSIVE ESTIMATES;OBS WIND 

1VAR=R[T] = °;2 VAR= + ;BIN ON R[T] 
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LOG AVERAGE PRED MSE PER BIN 



550 M3 V WIND;RECURSiVE ESTIMATES;OBS WIND 

1 VAR=R[T] = o;2 VAR= + ;BIN ON R[T] 
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LOG AVERAGE PRED MSE PER BIN 



250 MB U WINDiRECURSIVE ESTIMATES: 1 ST GUESS 
1VAR=WS[T]=o;2 VAR= + ;BIN ON WS[T] 
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250 MB V WIND;RECURSIVE ESTIMATES: 1 ST GUESS 
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Figure 6 



250 MB U WIND;RECURSIVE ESTIMATES; 1 ST GUESS 

1 VAR = R[T] = o;2 VAR = + ;BIN ON R[T] 
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LOG AVERAGE PRED MSE PER BIN 



250 MB V WINDiRECURSIVE ESTIMATES; 1 ST GUESS 

1 VAR=R[T] = o;2 VAR=+;BIN ON R[T] 
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LOG AVERAGE PRED MSE PER BIN 



500 MB U WiND;RECURS(VE ESTIMATES; 1 ST GUESS 

1 VAR=R[T] = o;2 VAR=+;BIN ON R[T] 
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LOG AVERAGE PRED MSE PER BIN 



500 MB V WiND;RECURSIVE ESTIMATES; 1 ST GUESS 

1 VAR=R[T] = o;2 VAR= + ;BIN ON R[T] 
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LOG AVERAGE PRED MSE PER BIN 
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